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Abstract 

Two  parameter  estimation  methods  are  presented  for  online  determination  of  parameter  values  using  a  simple  charge/discharge  model 
of  a  Sony  18650  lithium  ion  battery.  Loss  of  capacity  and  resistance  increase  are  both  included  in  the  model.  The  first  method  is  a  hybrid 
combination  of  batch  data  reconciliation  and  moving-horizon  parameter  estimation.  A  discussion  on  the  selection  of  tuning  parameters  for  this 
method  based  on  confidence  intervals  is  included.  The  second  method  uses  batch  data  reconciliation  followed  by  application  of  discrete  filtering 
of  the  resulting  parameters.  These  methods  are  demonstrated  using  cycling  data  from  an  experimental  cell  with  over  1600  charge-discharge 
cycles. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction  and  motivation 

The  overall  performance  of  batteries  deteriorates  over  time 
as  the  system  is  cycled  repeatedly  through  multiple  iterations 
of  charge  and  discharge.  For  high-performance  applications 
it  is  useful  to  have  accurate  knowledge  of  the  present  condi¬ 
tion  of  the  battery  as  well  as  the  remaining  battery  life.  The 
reduction  in  battery  performance  can  be  assumed  to  mani¬ 
fest  in  two  ways:  capacity  fade  display  in  the  reduction  in 
ability  of  the  battery  to  store  charge  and  increased  area  spe¬ 
cific  impedance  (ASI)  the  resistance  to  charge  transfer  which 
reduces  cell  potential. 

Given  the  importance  of  this  information  for  high  perfor¬ 
mance  electrochemical  systems,  significant  effort  has  been 
devoted  to  the  development  of  models  that  describe  the  dis¬ 
charge  behavior  of  batteries.  The  majority  of  these  models 
are  empirical  or  semi-empirical  at  best  [1,2],  having  limited 
predictive  capability  given  the  strong  dependence  of  cell  be¬ 
havior  on  factors  such  as  temperature  and  charge/discharge 
cycling  protocol.  However,  even  models  developed  using 
a  more  mechanistic  approach  [3],  are  susceptible  to  error. 
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These  models  are  typically  developed  using  data  that  is 
collected  under  tightly  controlled  experimental  conditions, 
which  simply  do  not  exist  in  real-world  applications.  Nev¬ 
ertheless,  these  models  can  be  regressed  to  data  and  can 
provide  an  accurate  representation  of  a  battery  for  a  few 
charging/discharging  cycles.  As  is  the  case  for  any  model 
prediction,  unknown  disturbances  and  unmodeled  phenom¬ 
ena  will  inevitably  influence  the  system,  leading  to  the  di¬ 
vergence  of  the  model  prediction  from  the  actual  battery 
performance. 

This  universal  drawback  of  modeling  is  typically  over¬ 
come  in  real-world  applications  (e.g.  weather  forecasting) 
by  periodically  updating  the  parameters  of  the  model  to  re¬ 
flect  the  current  state  of  the  system.  For  a  number  of  dy¬ 
namic  systems,  this  is  not  a  trivial  matter.  Often  it  is  dif¬ 
ficult  to  measure  these  parameter  values  without  disturb¬ 
ing  the  system,  if  it  is  even  possible  to  measure  them  di¬ 
rectly  at  all.  For  Li-ion  cells,  the  only  method  currently 
available  to  directly  measure  the  actual  capacity  is  to  stop 
the  cycling  and  physically  open  the  battery  in  an  inert  en¬ 
vironment.  Not  only  is  this  counterproductive  as  the  cy¬ 
cling  is  interrupted,  but  it  is  impossible  for  remote  appli¬ 
cations  such  as  satellite  power  systems.  Therefore,  there  is 
great  value  in  being  able  to  determine  updated  model  pa- 
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Nomenclature 

an  ith  element  on  the  diagonal  of  A 
an  nth  auto-regressive  coefficient  in  a  digital  filter 
A  parameter  covariance  matrix  in  parameter  es¬ 

timation 

bn  nth  moving-average  coefficient  in  a  digital  fil¬ 
ter 

bi  estimated  value  of  unknown  parameter  Pi 

dS®  arbitrary  process  disturbance  at  time/cycle  k 
e(k)  error  between  measured  and  model-predicted 
arbitrary  process  at  time/cycle  k\ 
i  current  drawn  from  battery  during  discharge 

(A) 

J  n  x  m  Jacobian  matrix  used  in  parameter  es¬ 

timation  for  the  rates  of  change  of  the  n  mea¬ 
surements  with  respect  to  the  m  undetermined 
parameters 

k  discrete-time  increment  or  alternatively  cycle 

number 

f  number  of  previous  cycles  considered  in  hybrid 

estimation  algorithm 

m  number  of  unknown  parameters  in  a  general 
parameter  estimation 

n  number  of  experimental  measurements  used  in 

a  general  parameter  estimation 
parameter  set  for  arbitrary  process  at 
time/cycle  k 

Q  battery  capacity  (A  h) 

R  internal  battery  resistance  (Q) 

s1  measurement  covariance  found  during  param¬ 

eter  estimation  that  approximates  true  covari¬ 
ance  a2 
t  time  (s) 

^cutoff  time  at  which  experimental  system  stops  dis¬ 
charging  based  on  a  cutoff  voltage  (s) 
t(  i-a/2)  signifance  statistic  at  (1  —  a)  confidence  level 
found  from  Student’s  t  distribution  with  v  de¬ 
grees  of  freedom 

u(k)  arbitrary  process  input  at  time/cycle  k 
U  empirical  correlation  between  SOC  and  open 
circuit  potential  of  the  battery  (V) 

V  experimentally  measured  battery  potential 
(V) 

V  model-predicted  battery  potential  (V) 

y(^)  arbitrary  process  measurement  at  time/cycle  k 
y(*)  model-predicted  value  of  arbitrary  process  at 

time/cycle  k 

z  unit  shift  operator 

Greek  symbols 

Pi  general  unknown  model  parameter  to  be  found 
by  estimation 


r Q 

move-weighting  penalty  on  prior  values  of  Q 
in  hybrid  estimation  algorithm 

rR 

move-weighting  penalty  on  prior  values  of  R 
in  hybrid  estimation  algorithm 

e 

state-of-charge  (SOC) 

^cutoff 

SOC  of  experimental  system  when  discharging 
stops  based  on  a  cutoff  voltage 

V 

number  of  degrees  of  freedom  in  a  parameter 
estimation,  calculated  by  n  —  m 

a2 

true  process  measurement  covariance  that  is 
approximated  by  s2  in  parameter  estimation 

0 

objective  function  used  in  hybrid  estimation  al¬ 
gorithm 

CJOc 

cutoff  frequency  for  digital  low-pass  filter  de¬ 
sign  (Hz) 

rameters  in  a  minimally  invasive  manner,  and  preferably 
online. 

In  the  following  work,  a  new  estimation  algorithm  —  com¬ 
bining  elements  of  both  batch  estimation  and  online  moving- 
horizon  estimation  [4-7]  —  is  proposed  to  keep  the  model 
parameters  updated  to  the  current  operating  state  of  the  bat¬ 
tery  on  a  cycle-to-cycle  basis.  Based  on  the  most  recent 
discharge  curve  of  a  cell,  capacity  and  resistance  are  cal¬ 
culated  through  the  minimization  of  an  error  function  that 
depends  on  the  mismatch  between  the  model  and  the  data 
from  the  present  cycle  as  well  as  a  weighted  penalty  for  the 
deviation  of  the  present  parameters  from  prior  values.  The 
addition  of  the  deviation  terms  smoothes  the  apparent  pa¬ 
rameter  drift  of  the  system  and  the  weightings  are  chosen 
to  balance  between  the  smoothness  of  the  parameter  drift 
and  the  achievement  of  the  smallest  possible  model/battery 
data  mismatch  for  a  given  cycle.  The  proposed  algorithm  is 
demonstrated  on  1600+  cycles  of  data  from  a  Sony  18650 
Li-ion  cell.  The  results  presented  here  are  presented  based 
on  several  different  tunings  of  the  parameters  of  the  algo¬ 
rithm  that  should  not  to  be  confused  with  the  model  param¬ 
eters.  As  an  alternative  method,  the  model  parameters  are 
obtained  using  successive  batch  estimations  for  each  cycle 
and  then  filtered  using  a  discrete,  recursive  low-pass  filter. 
The  results  from  the  filtering  analysis  are  compared  to  those 
from  the  new  algorithm,  and  analogies  are  drawn  between 
them. 

Regardless  of  the  method  chosen  to  analyze  the  data, 
it  is  apparent  that  the  parameters  do  not  obey  simple  dy¬ 
namics  that  can  be  predicted  a  priori.  Despite  being  the 
same  type  of  cells  running  the  same  protocol  simultane¬ 
ously  under  identical  conditions,  the  dynamics  of  the  ca¬ 
pacity  fade  and  increased  resistance  differ  significantly  be¬ 
tween  the  two  individual  experiments,  further  motivating  the 
need  for  this  online,  soft-sensing,  parameter  estimation  algo¬ 
rithm. 
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2.  Methodology 

The  proposed  algorithm  is  devised  as  a  method  to  observe 
slow  drifts  in  characteristics  of  Li-ion  batteries.  While  batch 
estimations  could  be  performed  for  each  cycle,  the  values  of 
the  parameters  obtained  in  this  manner  are  liable  to  fluctuate 
rapidly  in  a  relatively  uncorrelated  manner.  However,  it  is  ex¬ 
pected  that  these  model  characteristics  should  vary  smoothly 
and  gradually,  barring  the  occurrence  of  a  fault  in  the  sys¬ 
tem.  Therefore,  the  batch  estimations  are  performed  using  an 
objective  function  that  incorporates  parameter  values  from 
the  previous  cycle(s)  to  correlate  and  smooth  the  drift  of  the 
parameters. 

2.7.  Experimental  system 

Two  Sony  18650  Li-ion  batteries  were  used  to  collect  data 
over  multiple  charge-discharge  cycles.  These  batteries  were 
rated  for  1.4  Ah.  The  charging  protocol  was  as  follows:  con¬ 
stant  current  was  applied  at  1  A  until  the  battery  reached 
4.2  V,  then  constant  voltage  was  applied  at  4.2  V  until  cur¬ 
rent  dropped  to  50  mA.  The  discharge  protocol  was  as  fol¬ 
lows:  batteries  were  discharged  at  a  constant  1  C  rate  (1.4  A) 
then  stopped  when  the  voltage  reached  2.8  V.  The  discharged 
battery  then  waited  for  the  charger  to  become  available,  typ¬ 
ically  allowing  the  cell  to  rest  at  open  circuit  for  5-30  min. 
While  results  from  both  batteries  were  very  similar,  only  the 
results  from  the  first  battery  were  included  in  the  final  version 
of  this  work. 

A  schematic  of  the  data  collection  system  is  provided  in 
Fig.  1.  Charging  was  accomplished  using  a  Hewlett  Packard 
Model  663 2B  DC  Power  Supply.  An  Agilent  Technologies 
electronic  load,  Model  6060B,  was  used  to  draw  current  from 
the  batteries  on  discharge.  A  National  Instruments  General 
Purpose  Interface  Bus  (GPIB)  was  used  to  control  remotely 
the  load  and  power  supply.  Data  acquisition  and  control  were 
accomplished  using  Lab  VIEW®  software  version  6.1. 


2.2.  Dynamic  battery  discharge  model 

The  Matlab®  and  Simulink®  simulation  environments 
were  used  to  implement  this  estimation  routine.  A  model  of 
the  process  is  implemented  in  Simulink®.  A  simple  two- 
parameter  model  of  the  battery  is  defined  by  the  following 
equations: 


9 

U{6)  =  Yjapj  (2) 

j= 0 

v(t)  =  mm)  -  m ) 

where  0(t)  is  a  quantity  called  state-of-charge  (SOC)  that 
varies  between  0  and  1,  representing  fully  discharged  and 
fully  charged,  respectively.  The  time-dependent  behavior  of 
0(t)  is  given  by  Eq.  (1).  Battery  capacity  Q  is  assumed  to 
be  a  constant  parameter  during  the  course  of  a  single  bat¬ 
tery  discharge.  For  the  data  sets  used  in  this  study,  i(t),  or 
current  load  on  the  battery,  was  held  constant.  Thus,  6(f)  de¬ 
creases  linearly  over  time,  but  may  decease  at  different  rates 
as  the  model  capacity  Q  changes.  An  empirically  fit,  ninth- 
order  polynomial,  U ( 6 ),  maps  the  SOC  to  a  voltage  potential 
discharge  curve  being  produced  by  the  battery  (see  Eq.  (2)). 
However,  the  battery  itself  has  a  resistive  load,  R ,  so  the  ex¬ 
ternally  available  voltage,  V,  is  reduced  by  the  quantity  Rift) 
as  shown  in  Eq.  (3). 

The  model  parameters  from  Eq.  (2)  are  determined  from 
a  normal  discharge  of  the  battery  system  under  a  variety  of 
assumptions.  It  is  assumed  that  the  initial  capacity  of  the 
battery  is  as  stated,  1.4  Ah.  At  the  1  C  discharge  rate  (1.4  A) 
the  system  never  reaches  an  SOC  of  0,  since  discharge  of 
the  system  is  stopped  at  a  cutoff  voltage  of  2.8  V.  Given  a 


Electric  Load 


Mechanical  Relays 


Fig.  1.  Schematic  layout  of  battery  charge  discharge  system. 
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constant  discharge  rate,  i,  the  SOC  at  cutoff  is  given  by 


^cutoff  —  1 


^cutoff 


(4) 


Assuming  the  initial  system  resistance  R  is  0.3  £2,  the  result¬ 
ing  voltage  profile  V(t )  for  a  single  discharge  at  a  1  C  rate 
can  then  be  fit  to  a  high-order  polynomial.  Using  the  first 
discharge  cycle  for  the  battery,  the  following  relationship  is 
found: 


2.3.1.  Hybrid  estimation  implementation 

The  online  estimation  routine  consisted  of  a  nonlinear  op¬ 
timization  at  each  analyzed  cycle  solved  by  the  unconstrained 
optimization  algorithm  fminsearch  (Nelder-Mead)  provided 
by  Matlab®.  The  objective  function  used  was: 

Uf(y-V)2  fa  i 

<P(G.-.  Ri)  =  —t — - +  rQ  2~j(Qi  -  Qi-j f 

tf-to  u 


U{0)  =  4  -  2510  +  266 02  -  1352<93  +  41486>4  -  8.O7305 
+  994606  -  7472<9731O708  -  54309  (5) 

Note  that  this  model  is  not  based  on  the  open  circuit  potential. 
Also,  the  polynomial  is  never  updated  over  the  course  of 
the  estimation  procedure.  The  values  for  Q  and  R  will  be 
determined  in  the  estimation  procedures. 

2.3.  Hybrid  estimation  methodology 

The  proposed  estimation  algorithm  can  be  viewed  as  an 
extension  of  traditional  online  parameter  estimation.  First, 
the  known  input  values  for  the  Mi  cycle,  uSk\  the  current,  are 
used  to  force  the  battery  to  charge  and  discharge.  Unmea¬ 
sured  disturbances,  d^k\  influence  the  output  of  the  battery  in 
cycle  k.  These  disturbances  could  include  temperature  vari¬ 
ation,  unmodeled  dynamics,  or  fault  occurrences  in  the  bat¬ 
tery.  Disturbances  affect  the  measured  voltage  profile  for  the 
discharge  cycle,  y^k\  but  do  not  affect  the  modeled  output, 
y(k\  The  model  discharge  profile  also  depends  on  the  current 
model  parameters,  P^k\  where  =  [Q^  R^]T  .  The  er¬ 
ror  between  the  battery  and  the  model,  e^k\  is  then  passed 
to  the  optimization  routine  to  determine  the  best  parameter 
values  for  the  current  cycle. 

The  optimization  routine  attempts  to  minimize  a  weighted 
combination  of  the  mean- squared  error  (MSE)  for  the  battery 
and  a  penalty  based  on  the  squared  value  of  deviations  from 
previous  parameter  values.  For  this  application,  MSE  is  cho¬ 
sen  as  the  model  —  data  error  metric  over  a  simpler  form  such 
as  sum-squared  error  (SSE).  Since  the  length  of  a  discharge 
cycle  is  not  constant,  but  rather  is  determined  when  the  sys¬ 
tem  reaches  the  cutoff  voltage,  the  number  of  data  points  in 
each  cycle  is  variable.  Hence,  the  value  of  an  error  metric 
such  as  SSE  is  not  directly  comparable  between  cycles  of 
different  length.  However,  the  averaging  process  involved  in 
MSE  eliminates  the  discrepancies  between  cycles  of  different 
length,  making  it  a  more  appropriate  metric  for  this  system. 

The  minimization  of  the  objective  function  is  accom¬ 
plished  by  adjusting  the  values  of  P^k\  Each  cost  function 
evaluation  of  the  optimization  returns  a  new  set  of  values  for 
to  the  model,  which  then  recalculates  y^  for  the  cycle. 
Once  the  optimization  terminates,  the  final  values  for  P^are 
output  and  also  saved  for  use  in  subsequent  estimation  itera¬ 
tions  for  later  cycles.  The  algorithm  increments  itself  to  the 
next  cycle  number  and  repeats  the  entire  process. 


+  rRJ22~j(Ri-  Ri-j)2  (6) 

7=0 

where  Qi  and  Rf  are  the  model  parameters  being  fit  for  cycle 
i,  and  V  and  V  are  the  battery  modeled  and  experimental  volt¬ 
ages,  respectively.  The  initial  and  final  times  of  the  data  set 
are  to  and  tf.  Since  the  data  is  sampled  at  discrete  intervals, 
the  data  is  linearly  interpolated  and  the  continuous  value  of 
the  model  —  data  error  is  integrated  by  Simulink.  Qi-j  and 
Ri-j  are  the  values  of  Q  and  R  from  the  (i  —  j) th  cycle,  and 
Fq  and  Ur  are  the  relative  weightings  of  the  move-penalty 
terms  of  the  objective  function.  Fq  and  Fr  are  the  magni¬ 
tudes  of  the  weights  for  the  move  penalty  for  the  (/  —  l)th 
cycle,  and  the  weighting  for  each  cycle  prior  to  that  is  re¬ 
duced  by  a  factor  of  2.  There  are  two  components  to  this 
objective  function:  the  integral  term  measuring  the  ability  of 
the  model  to  reproduce  the  data  and  the  terms  that  penalize 
large  deviations  in  parameter  values  between  runs.  It  is  nec¬ 
essary  to  choose  the  Fq  and  Fr  terms  carefully  so  that  they 
do  not  dominate  the  model  error  term  in  Eq.  (6).  The  ratio 
between  Fq  and  Fr  is  also  important,  allowing  both  param¬ 
eters  to  have  similar  relative  impact  even  though  Q  and  R  are 
of  different  magnitudes. 

In  addition  to  the  weighting  parameters  Fq  and  Fr,  an¬ 
other  adjustable  parameter  for  tuning  algorithm  performance 
is  i,  the  size  of  the  weighting  horizon  in  number  of  cycles. 
The  case  where  t  equals  0  is  a  special  case,  since  the  move- 
penalty  terms  of  the  objective  function  become  0,  and  the 
objective  function  simply  becomes  an  unconstrained  mean- 
squared  error  (MSE)  for  the  present  discharge  cycle.  Due  to 
the  geometrically  decreasing  weights  on  move  penalties  at 
larger  values  of  i,  the  incremental  effects  of  increasing  i  to 
values  larger  than  4  are  essentially  negligible. 

2.3.2.  Parameter  confidence  intervals  from  hybrid 
estimation 

In  many  cases,  it  is  useful  to  obtain  confidence  intervals 
for  regressed  parameters  as  a  measure  of  the  reliability  of 
the  results.  A  straightforward  statistical  technique  exists  for 
computing  these  intervals  for  simple  linear  least  squares,  as¬ 
suming  that  all  errors  are  normally  distributed.  A  procedure  is 
also  available  for  unconstrained  nonlinear  least  squares  (un¬ 
der  the  same  error  assumptions),  but  it  is  significantly  more 
involved  than  linear  least  squares.  Given  the  unconventional 
objective  function  used  in  the  hybrid  estimation,  it  would  ap¬ 
pear  difficult  to  derive  such  a  procedure  for  this  new  method. 
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However,  by  slightly  reformulating  the  model-data  error  cal¬ 
culation  and  considering  the  parameter  drifts  as  additional 
measurements,  it  will  be  shown  that  a  method  exits  that  can 
be  used  to  compute  confidence  intervals. 

Constantinides  and  Mostouh  [8]  discuss  the  calculation  of 
confidence  intervals  from  nonlinear  least  squares  regressions 
in  which  more  than  one  type  of  data  is  used.  In  this  case, 
parameters  are  chosen  to  minimize  a  weighted  SSE  objec¬ 
tive  function.  For  the  general  case  with  v  different  types  of 
measurements  Yj,  the  objective  function  to  be  minimized  is 
given  by 

V 

&  =  Ewj(Yj-yjA(Yj-yj)  (7) 

j=  1 

In  the  case  where  the  variance  for  each  type  of  measure¬ 
ment  crj  is  known,  then  the  weighting  factors  wj  can  be  com¬ 
puted  directly  and  are  inversely  proportional  to  crj.  However, 
it  is  not  uncommon  for  the  variances  to  be  unknown,  so  one 
must  estimate  values. 

Now,  consider  the  case  of  the  battery  system.  In  order  to 
use  this  approach,  the  error  between  the  battery  and  model 
will  need  to  be  handled  as  SSE  instead  of  MSE,  despite  the 
previously  mentioned  advantages  of  MSE.  It  also  appears 
that  there  is  only  one  measurement:  the  voltage  in  the  bat¬ 
tery.  However,  the  move-penalty  terms  that  are  included  in 
the  original  objective  function  of  Eq.  (6)  can  be  considered 
additional  measurements.  In  the  case  where  1  =  1,  there  end 
up  being  three  measurements: 


Yi  =  V,  Y2  =  (Qi  -  Qi-\),  Y3  =  -  Ri-i) 


(8) 


with  a  similar  expression  for  the  /?/  term.  If  this  is  truly  valid, 
then  the  objective  function  further  simplifies  to 

n 

=  Y,{%  ~  +  re(2i  -  Qi- 1)2  +  Wi  -  i)2 

k=l 

(12) 

which  is  nearly  identical  to  the  objective  function  proposed 
in  this  work  for  the  case  where  l  =  1 .  In  a  future  case  where 
a  model  or  proposed  model  of  the  time-dependent  drift  of  Q 
and  R  is  known  (such  as  a  square  root  dependence  on  time), 
then  that  model  could  alternately  be  used  to  find  values  for 
(Qi  ~  Qi- 1)  and  (Ri  —  Ri-Q.  Confidence  intervals  for  esti¬ 
mated  parameters  are  given  by  Constantinides  and  Mostouh 
in  the  form 


bi  )S\f&U  —  Pi  ^  bi  H-  t(\—a/2)S\f(2ii  (13) 


where  bi  is  the  value  of  the  parameter  Pi  obtained  by  mini¬ 
mizing  0,  t( i-a/2)  the  significance  statistic  for  a  confidence 
level  of  (1  —  ot\  s  the  square  root  of  the  approximate  vari¬ 
ance  of  the  error  and  an  the  ith  element  of  the  diagonal  of  the 
parameter  covariance  matrix.  Values  of  t( \-a/2)  can  either  be 
found  by  lookup  in  statistical  tables  of  the  Student’s  t  distri¬ 
bution  or  direct  computation  for  a  given  number  of  degrees 
of  freedom  v  and  the  desired  significance  level,  v  is  given  by 
the  number  of  measurements  less  the  number  of  parameters 
being  fit.  When  v  becomes  large,  this  statistic  approaches  that 
of  the  normal  distribution.  Values  for  s  are  obtained  from  the 
minimum  value  of  the  objective  function: 


v 


(14) 


Thus,  the  objective  function  becomes 


Finally, 


<*>,■ = yy  v*  -  vo + rQ[(Qf^Qi- 1)  -  (Qt  -  i)]2 

k=  i 


+  I)  -  (Ri  -  Ri-i)]2  (9) 


A  = 


j= 1 


where 


where 

w\  =  1,  w2  =  rQ ,  w3  =  rR 


(10) 


9^ i  ...  9>7.i 
db\  dbm 


(15) 


noting  that  w\  is  assumed  to  be  1  and  the  arbitrary  weighting 
parameters  Tq  and  TR  are  chosen,  because  nothing  is  known 
a  priori  about  the  variance  of  the  two  “pseudo-measurements” 
relative  to  the  measurement  error.  Since  no  assumptions  are 
made  as  to  how  the  parameters  Q  and  R  drift  with  time,  one 
hypothesis  is  to  consider  that  Q  and  R  are  Gaussian  random 
variables  and  that  the  values  of  Qi  and  Ri  found  in  each  cycle 
are  samples  from  those  distributions.  In  this  case,  the  model 
of  the  cycle-to-cycle  drift  of  Q  is  given  by 


(ID 


Jj  = 


d_hn 

.  dbi 


d_hn 
dbm  _ 


(16) 


for  the  case  where  there  are  n  instances  of  the  jth  measure¬ 
ment  and  m  parameters  being  fit.  Now  consider  the  case  of 
battery  discharge.  The  calculation  of  s  is  very  simple.  Using 
the  simplified  objective  function  for  the  problem,  one  obtains 


2  <Pi(Qi,Ri) 

s  =  - 


<Pi(Qi,Ri) 

n  T  2  —  2 


Ri) 


(Qi  -  Qi- 1)  =  E[Qi  -  Qi- 1]  =  0 


v 


n 


(17) 
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where  again,  n  is  the  number  of  voltage  measurements.  Mea¬ 
surements  Y2  and  Y3  have  very  simple  Jacobians  to  compute: 


jJj2  =  [1  0]T[1  0] 


1  0 
0  0 


(18) 


Hh  =  [0  i]'[0  i]  = 


0  0 
0  1 


(19) 


The  Jacobian  for  the  voltage  measurements  is  somewhat  more 
involved.  The  voltage  model  for  constant-load  discharge  is 
given  by  Eqs.  (l)-(3).  Therefore,  the  Jacobian  of  the  voltage 
is  given  by 


dV 

do 

dV 

~d0 

0=0(ti)  9 Q 

,=n  ™ 

J 1  = 


dv 

de 

dv 

~de 

0=6(t„)  dQ 

t=tn  ™ 

which  simplifies  somewhat  to 
itj 


Ji  = 


dV 

~d0 


dv 

Ho 


0=0(t  1) 


e2 


e=e(t„) 


l  tn 

G2 


—  l 


—  l 


(20) 


(21) 


Using  the  definitions  for  J  the  value  for  A  is  obtained: 


A  =  [Jj  JX  +  rQf  h  +  rRf  J3]  1 


t  (S)2 

k=l 

n 

Emm 

dQ  dR 


Emm 

dQ  dR 


i:(§)2 


-1 


0 

rR 


(22) 


Using  these  equations,  it  is  now  possible  to  compute  the 
confidence  intervals.  Note  that  Tq  and  Tr  both  influence  the 
confidence  intervals  by  their  inclusion  in  A.  Since  Tq  and  7^ 
are  not  based  on  actual  variances,  it  must  be  emphasized  that 
excessively  large  values  of  these  parameters  will  yield  nar¬ 
rower  confidence  intervals  than  what  may  actually  be  valid.  In 
practice,  this  procedure  calculates  confidence  intervals  very 
similar  to  those  that  would  be  obtained  if  the  move-penalty 
terms  were  not  incorporated.  As  seen  in  Fig.  4,  the  confidence 
range  begins  decreasing  with  increasing  Tq  and  Tr.  Since 
larger  values  of  Tq  and  7#  correspond  to  smaller  variances  or 
expected  drifts  in  Q  and  R ,  it  is  expected  that  one  would  have 
a  higher  degree  of  confidence  in  parameters  that  are  known 
to  drift  less.  Nevertheless,  values  of  Tq  and  Tr  should  be 
chosen  judiciously. 


2.4.  Digital  filtering 

A  second  analysis  method  is  proposed  where  the  model 
parameters  are  calculated  at  each  cycle  using  simple  batch 
estimation  techniques.  Then  a  discrete,  low-pass  filter  is  ap¬ 
plied  to  the  values  to  remove  the  high-frequency  variations. 
Due  to  its  simplicity,  a  second-order  Butterworth  filter  (HR) 
was  constructed  using  the  Matlab®  Signal  Processing  Tool¬ 
box.  When  developing  a  low-pass  filter,  a  crossover  frequency 
must  be  supplied.  For  a  discrete  (digital)  filter,  it  is  gener¬ 
ally  sufficient  to  specify  the  ratio  between  the  cutoff  fre¬ 
quency  and  the  sampling  rate.  The  battery  system  is  not  truly 
a  discrete-time  system,  since  the  data  points,  i.e.  estimated 
parameter  values,  do  not  correspond  to  sequential,  evenly 
spaced  samples.  Instead,  each  estimate  corresponds  to  an  en¬ 
tire  discharge  cycle,  which  does  not  have  a  fixed  duration. 
However,  for  the  purposes  of  this  analysis,  the  cycle-to-cycle 
frequency  is  arbitrarily  assumed  to  be  1  Hz.  Therefore,  a  fil¬ 
ter  designed  with  a  cutoff  frequency  of  0.01  Hz  will  reject 
dynamics  of  parameter  changes  that  occur  over  fewer  than 
100  cycles.  In  other  words,  the  filtered  values  for  capacity 
and  resistance  should  only  exhibit  dynamics  that  occur  over 
at  least  100  cycles. 


3.  Results 

Several  different  studies  were  performed  to  gauge  the  per¬ 
formance  of  the  proposed  hybrid  estimation  algorithm.  One 
of  the  main  assumptions  in  this  problem  is  that  the  param¬ 
eters  drift  slowly  enough  that  they  can  be  held  constant  for 
the  duration  of  a  single  cycle.  If  this  assumption  is  valid,  it  is 
improbable  that  the  actual  parameter  value  will  have  changed 
significantly  from  one  cycle  to  the  next.  Therefore,  it  is  de¬ 
sirable  to  tune  the  algorithm  to  produce  parameter  profiles 
with  reduced  high-frequency  “chatter”  in  the  signal,  much 
like  a  low-pass  filter.  However,  if  the  move-penalty  weights 
are  too  high,  the  parameter  values  will  lag  and  the  error  in 
the  resulting  model  could  be  unacceptable.  The  studies  dis¬ 
cussed  herein  consider  the  effects  of  changing  move-penalty 
weights  and  varying  the  length  of  the  move-penalty  horizon. 
These  results  are  compared  to  results  obtained  by  digital  fil¬ 
tering.  From  these  studies,  recommendations  are  made  for 
appropriate  magnitudes  of  algorithm  tuning  parameters  and 
filter  specifications. 

3.1.  Effects  of  hybrid  estimation  algorithm  tuning 
parameters 

One  experiment  was  performed  where  the  move-penalty 
weights  r<2  and  Fr  were  varied  while  the  move  horizon 
length  i  was  fixed  at  2.  The  results  of  this  experiment  are 
presented  in  Fig.  2.  As  was  expected,  the  smallest  weights 
(Eg  =  2  and  7#  =  4)  produced  results  that  were  extremely 
similar  to  the  unweighted  case.  The  estimated  values  of  Q 
and  R  exhibit  a  very  high  variation  over  an  underlying  steady 
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Fig.  2.  Capacity  fade  and  resistance  change  with  various  values  of  Cg  and  7^.  Note  that  =  2Tq.  The  move  horizon  length  parameter  i  was  fixed  at  2. 


decrease  and  increase,  respectively.  It  is  doubtful  that  the  ac¬ 
tual  condition  of  the  battery  changes  as  rapidly  as  predicted. 
Alternately,  the  case  with  the  largest  weights  ( Tq  =50  and 
rR  =  100)  yields  parameter  trends  that  are  very  smooth 
and  gradual,  but  appear  to  lag  significantly  behind  the  un¬ 
weighted  results.  This  too  is  unacceptable.  Finally,  a  reason¬ 
able  compromise  weighting  appears  to  be  about  ( Fg  =  10 
and  rR  =  20).  In  this  case,  the  highest  frequency  noise  com¬ 


ponents  are  not  present,  but  the  estimated  values  still  converge 
to  the  long-term  trends  of  the  unweighted  case  rapidly. 

A  second  study  was  conducted  varying  the  move  hori¬ 
zon  length,  i,  while  maintaining  the  move-penalty  weights 
at  the  values  determined  previously  to  be  most  reasonable 
(rg  =  10  and  rR  =  20).  As  seen  in  Fig.  3,  the  effect  of  the 
increasing  horizon  length  was  to  smooth  out  the  predicted 
value  curve  at  the  expense  of  convergence  rate,  similar  to 
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Fig.  3.  Capacity  fade  and  resistance  change  using  various  values  of  t  in  the  estimation.  The  move-penalty  weighting  parameters  were  fixed  at  Tq  =  10  and 
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Fig.  4.  95%  Confidence  intervals  for  Q  and  R  over  1600  charge-discharge 
cycles. 


increasing  the  weighting  factors  Fq  and  Fr.  Fig.  3  also  illus¬ 
trates  that  the  differential  effect  between  1  —  4  and  i  =  2  is 
approximately  equivalent  to  that  between  i  =  2  and  i  =  1 . 
With  the  weighting  scheme  currently  implemented,  the  ef¬ 
fects  of  setting  i  larger  than  4  are  essentially  negligible.  For 
the  particular  values  of  the  weighting  parameters  used  in  this 
study,  a  value  of  i  =  2  yields  the  best  blend  of  smoothing  the 
profiles  without  adding  excess  lag  and  unduly  reducing  the 
accuracy  of  the  model. 

Fig.  4 

Finally,  there  is  no  utility  in  adjusting  the  algorithm  pa¬ 
rameters  Fq,  Fr ,  and  i  to  obtain  smooth  parameter  trends 
if  it  produces  unacceptable  error  between  the  model  and  the 
battery  voltage.  Therefore,  the  hybrid  estimation  was  per¬ 
formed  using  the  ideal  parameter  values  determined  from  the 
previous  two  experiments.  Three  discharge  cycles  were  se¬ 
lected  from  early,  mid,  and  late  in  the  data:  cycles  1,  800,  and 
1425,  respectively.  The  comparison  between  the  data  and  the 
model  using  the  estimated  parameters  is  presented  in  Fig.  5. 
Clearly,  the  agreement  is  quite  good,  particularly  for  cycles  1 
and  800.  Although  the  results  for  cycle  1425  were  not  quite 
as  accurate,  they  are  still  encouraging,  considering  the  sim¬ 
plistic  empirical  nature  of  the  underlying  model. 

3.2.  Digital  filtering 

Filtering  has  been  suggested  as  an  alternate  means  of  uti¬ 
lizing  noisy  parameter  estimates  for  individual  cycles  and 
producing  smooth  trajectories  for  parameter  values  over  time. 
In  practice,  it  has  proved  to  be  a  very  simple  and  effec¬ 
tive  method  for  accomplishing  this  task.  Using  readily  avail¬ 
able  tools  such  as  Matlab®  Signal  Processing  Toolbox  and 
Simulink®,  a  number  of  Butterworth  low-pass  digital  filters 
were  synthesized  with  various  crossover  frequencies  and  used 
to  process  model  parameter  values  obtained  from  estimations 
on  single  cycles.  The  general  form  of  the  second-order  recur¬ 


Fig.  5.  Model  fidelity  for  the  battery  for  three  cycles  {1,  800,  1425}  using 
parameters  obtained  through  hybrid  estimation  and  low-pass  filtering. 


sive  digital  filter  (HR)  is  given  by 

y(k)  =  -a\y(k  -  1)  -  a2y(k  -  2)  +  b0y(k)  +  b\y(k  -  1) 

+  b2y(k-  2)  (23) 


where  y  is  the  filtered  value  of  the  quantity  y  and  k  is  the 
discrete-time  unit.  The  signs  on  the  coefficients  are  defined  in 
this  manner,  because  the  digital  filter  can  also  be  equivalently 
represented  in  the  following  discrete-time  transfer  function: 


y(z)  = 


bo  +b\z  1  +b2z  2 

t-, — =r~; — 

1  +  a\z  1  +  a2z  z 


(24) 


where  z  is  the  unit  shift  operator.  The  results  are  shown  in 
Fig.  6. 

In  Fig.  6,  the  cutoff  frequencies  used  for  the  various  fil¬ 
ters  were  {0.001,0.01,0.1}  Hz  for  the  capacities,  Q ,  and 
{0.0005,  0.005,  0.05}  Hz  for  the  resistances,  R.  The  filter  co¬ 
efficients  for  the  Q  values  are  listed  in  Table  1 ,  while  those  for 
R  are  listed  in  Table  2.  Based  on  the  behavior  of  these  filters, 
cutoff  frequencies  of  approximately  10-2  Hz  are  appropriate 
for  this  system. 

To  ensure  that  parameter  values  obtained  by  the  filters 
were  still  valid  for  the  system,  these  filtered  parameter  val¬ 
ues  were  used  to  generate  model  output  for  cycles  1,  800, 
and  1425,  as  was  done  for  the  hybrid  estimation  results.  The 
model  output  was  compared  to  the  experimental  output  in 


Table  1 

Filter  coefficients  (as  defined  in  Eq.  (23))  used  for  filtering  Q 


Coefficient 

coc  =  0.001 

a)c  =  0.01 

CDc  =0.1 

a\ 

—  1.991e  +  00 

-1.911c +  00 

-1.143c +  00 

a2 

9.912c -01 

9.150e  —  01 

4.128c -01 

bo 

9.826e  -  06 

9.441e  -  04 

6.745e  -  02 

bx 

1.965c -05 

1.889c -03 

1.349c  —  01 

b2 

9.826e  -  06 

9.441e  -  04 

6.745e  -  02 

coc  is  design  cutoff  frequency  in  Hz.  Values  are  obtained  using  the  Butter- 
worth  filter  tool  in  the  Matlab  Signal  Processing  Toolbox. 
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Fig.  6.  Capacity  fade  and  resistance  change  in  the  battery  smoothed  by  three  digital  Butterworth  filters. 


Table  2 

Filter  coefficients  (as  defined  in  Eq.  (23))  used  for  filtering  R 


Coefficient 

coc  =  0.0005 

coc  =  0.005 

coc  =  0.05 

a\ 

-1.9966* +  00 

—  1.9566? +  00 

-1.5616? +  00 

a2 

9.956e  -  01 

9.565e  -  01 

6.4136*  -  01 

bo 

2.462e  -  06 

2.4146*  -  04 

2.0086?  -  02 

bi 

4.924e  -  06 

4.8276*  -  04 

4.0166*  -  02 

bi 

2.462e  -  06 

2.4146*  -  04 

2.0086?  -  02 

co c  is  design  cutoff  frequency  in  Hz.  Values  are  obtained  using  the  Butter- 
worth  filter  tool  in  the  Matlab  Signal  Processing  Toolbox. 


Fig.  5.  There  appears  to  be  good  agreement  between  the 
model  using  filtered  parameters  and  experiment  over  all  cy¬ 
cles.  Moreover,  the  parameters  obtained  using  this  method 
appear  to  offer  comparable  performance  to  those  obtained 
using  the  hybrid  estimation  technique. 


4.  Discussion 

It  is  well  known  that  capacity  of  electrochemical  cells 
decreases  with  time  [1,3].  A  battery  that  has  been  recharged 
100  times  cannot  store  as  much  power  as  a  brand  new  battery. 
Consequently,  it  was  expected  that  the  estimation  algorithm 
would  show  the  capacity  parameter  in  the  model  studied  de¬ 
creasing  as  the  cycle  number  increased.  Additionally,  due  to 
irreversible  reactions,  the  internal  resistance  of  a  battery  will 
increase  with  the  number  of  cycles.  Both  of  these  experi¬ 
mentally  observed  trends  were  captured  by  the  tested  hybrid 
estimation  algorithm  and  discrete  filtering  of  batch  estimation 
results. 

The  new  algorithm  was  designed  as  a  hybrid  between  two 
more  traditional  types  of  estimation:  offline  batch  estima¬ 


tion  and  online  recursive  estimation.  Each  estimation  so¬ 
lution  found  is  a  result  of  the  solution  of  a  batch  estima¬ 
tion  problem.  However,  knowledge  gained  from  solutions 
of  the  problem  at  earlier  cycles  is  used  to  weight  the  ob¬ 
jective  function  used  in  the  current  estimation,  thereby  in¬ 
creasing  the  confidence  in  the  parameter  estimates,  more  akin 
to  recursive  estimation.  This  approach  assumes  that  the  dy¬ 
namics  of  the  parameters  being  estimated  are  slow  with  re¬ 
spect  to  the  length  of  each  data  cycle.  Hence,  they  are  as¬ 
sumed  constant  for  the  duration  of  a  single  cycle,  but  can 
vary  from  cycle-to-cycle.  See  [9]  for  a  discussion  of  cate¬ 
gories  of  adjustable  parameters.  The  hybridization  of  batch 
and  recursive  estimation  is  a  novel  approach  to  data  recon¬ 
ciliation. 

One  typical  downfall  of  batch  estimation  for  use  in  on¬ 
line  applications  is  that  the  algorithms  used  are  typically  not 
fast  enough  to  achieve  a  solution  within  one  time  interval. 
This  technique  mitigates  this  issue,  since  a  solution  is  com¬ 
puted  once  per  cycle  and  not  once  per  sample  period.  The 
battery  system  used  to  generate  data  for  this  study  has  a  sam¬ 
ple  time  of  1  s,  but  a  complete  cycle  is  about  50  min  long. 
The  solution  times  of  the  estimations  are  approximately  15  s, 
making  it  impossible  to  solve  this  type  of  problem  at  each 
time  step,  but  very  practical  on  a  cycle-to-cycle  basis.  This 
technique  should  be  applicable  to  a  number  of  batch  pro¬ 
cesses.  Any  process  that  operates  in  repeated,  discrete  cycles 
on  the  order  of  several  minutes  or  higher  may  be  able  to  ben¬ 
efit. 

The  theoretical  mathematical  properties  of  this  algorithm 
have  not  yet  been  analyzed.  Investigation  of  the  stability,  con¬ 
vergence,  robustness  in  the  presence  of  plant/measurement 
uncertainty,  or  any  other  desirable  property  would  be  a  non¬ 
trivial  exercise  at  best.  It  is  certainly  an  area  open  for  future 
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research.  This  algorithm  has  provided  excellent  performance 
in  the  application  problem  studied,  despite  the  absence  of 
rigorous  proofs  of  stability  or  convergence.  The  absence  of 
such  proofs  surrounding  an  algorithm  should  not  be  the  sole 
criterion  for  its  dismissal;  Allgower  et  al.  [7]  state  that  the 
extended  Kalman  filter  (EKF)  does  not  have  a  lot  of  strong 
theory  supporting  its  use,  but  has  nevertheless  proved  use¬ 
ful  in  a  number  of  applications.  The  EKF  has  also  been  used 
previously  for  online  estimation  of  battery  capacity  fade  [10- 
12].  However,  it  will  not  work  for  the  battery  model  used  in 
this  discussion,  since  this  model  is  not  observable  in  the  tra¬ 
ditional  dynamical  systems  context,  a  strict  requirement  of 
the  Kalman  filter  and  EKF.  Observability  for  linear  systems 
is  discussed  in  both  [13,14]  (among  others),  with  a  signif¬ 
icant  discussion  of  the  Kalman  filter  included  in  [13].  The 
concept  of  observability  is  extended  for  nonlinear  systems 
in  [15].  In  addition  to  observability,  the  extended  Kalman 
filter  requires  linearization  of  the  nonlinear  dynamic  model 
and  estimates  of  the  state  and  measurement  noise  characteris¬ 
tics.  Depending  on  the  system,  these  requirements  can  prove 
unwieldy.  Additionally,  the  EKF  and  other  online  recursive 
estimation  techniques  compute  updates  at  each  sample  point. 
For  a  system  such  as  the  battery  studied  here,  sampling  occurs 
frequently  (on  the  order  of  seconds),  whereas  the  parameter 
dynamics  of  interest  occur  on  a  much  longer  time  scale  (on 
the  order  of  days,  weeks,  or  months).  Consequently,  even  if 
it  were  possible  to  use  an  EKF,  it  would  require  substantially 
more  computational  effort  for  very  little  added  benefit.  Thus, 
the  hybrid  estimation  method  proposed  here  is  useful  and 
advantageous  in  this  circumstance  since  it  offers  the  “fading 
memory”  effects  of  recursive  estimation  at  a  more  reasonable 
intermediate  time  scale. 

Digital  filtering  also  appears  to  be  a  viable  approach  for 
developing  smooth  parameter  estimates  in  this  electrochem¬ 
ical  system.  Filter  tunings  are  readily  obtained  and  simple 
to  implement.  Using  the  filters,  it  is  possible  to  obtain  pa¬ 
rameter  profiles  that  are  relatively  smooth  without  exhibiting 
excessive  lag.  Therefore,  model  fidelity  remains  high.  Ad¬ 
ditionally,  the  filtering  work  provides  insight  on  methods  to 
further  improve  the  hybrid  estimation  algorithm.  Since  it  is 
assumed  that  the  process  and/or  measurement  are  subject  to 
significant  noise,  the  recursive  low-pass  filters  place  a  low 
weight  on  the  present  measured  value  relative  to  prior  filtered 
estimates.  This  weighting  scheme  suggests  that  the  current 
scheme  employed  for  the  weighting  parameters  through  the 
move  horizon  may  not  be  an  appropriate  one. 

Nevertheless,  both  methods  demonstrate  promising  means 
to  monitor  unmeasurable  parameters  of  electrochemical  sys¬ 
tems  online  without  relying  on  semi-empirical  correlations 
developed  under  highly  controlled  laboratory  conditions. 

5.  Conclusions 

Two  methods  have  been  presented  for  analyzing  capacity 
fade  in  Fi-ion  batteries.  Both  the  hybrid  estimation  and  the 


digital  filtering  techniques  offer  the  user  the  ability  to  moni¬ 
tor  the  long-term  condition  of  the  battery  while  minimizing 
the  noise  attributable  to  cycle-to-cycle  variation.  More  tradi¬ 
tional  approaches  for  online  parameter  tracking,  particularly 
the  extended  Kalman  filter,  are  impractical  for  this  type  of  ap¬ 
plication  given  the  orders  of  magnitude  difference  between 
sampling  interval  and  the  time  scales  of  parameter  drift.  Ad¬ 
ditionally,  the  EKF  is  precluded  for  this  particular  application 
since  the  dynamic  battery  model  used  is  not  observable.  Fi¬ 
nally,  results  from  the  hybrid  estimation  method  can  be  used 
under  certain  conditions  to  compute  confidence  intervals  for 
the  unknown  parameters  Q  and  R ,  if  a  quantitative  measure 
of  reliability  is  desired. 
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